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We present a theoretical study of the quantum depletion of microcavity polaritons that are excited 
with a resonant laser pulse. The dynamics of the quantum fluctuations are interpreted in the con¬ 
text of quantum quenches in general and in terms of the dynamical Casimir effect in particular. We 
compute the time evolution of the first and second order correlation functions of the polariton con¬ 
densate. Our theoretical modelling is based on the truncated Wigner approximation for interacting 
Bose gases. For homogeneous systems, analytical results are obtained in the linearised Bogoliubov 
approximation. Inhomogeneous systems are studied numerically by Monte Carlo simulations. 


I. INTRODUCTION 

Interaction quenches in quantum many body systems 
have become an active research fielcpJ^I, mainly thanks to 
the great degree of controllability of ultracold atoms with 
Feshbach resonances and optical latticed. Recently, a 
complementary platform for quantum many body physics 
that has been developed, namely exciton-polariton quan¬ 
tum fluid^. The distinctive characteristic of these sys¬ 
tems is that the polaritons are a superposition of light 
and matter excitations. A first advantage of the light 
component is that it allows for a straightforward diag¬ 
nostic of the fluid by means of standard quantum optical 
techniques. A second advantage is that polaritons can 
be created by an external laser field. It is this feature 
of polariton condensates that is of particular interest in 
the context of quantum quenches, since it allows to im¬ 
plement a sudden change in the many body system. 

The situation that we will consider here is an instanta¬ 
neous injection of polaritons in a coherent stat^. Since 
this is not the ground state of the interacting many-body 
system, a non-trivial time evolution will result. In our 
previous work, we have shown that a dynamical Casimir 
effec'P® takes place in terms of the Bogoliubov excita¬ 
tions on top of the coherently created polariton state. 
Indeed, the sudden creation of a condensate quenches 
the vacuum from the trivial one to the Bogoliubov vac¬ 
uum, resulting in an excitation of the system. Part of 
the motivation for the study of the dynamical Casimir 
effect stems from connections with the Hawking-Unruh 
effect, whose sonic versiorP^ is getting within the reach 
of experiments with polariton^ and ultracold atomJ^^. 

The analogy with an interaction quench in cold atom 
systems is direct, since our proposal is equivalent to a 
sudden increase of the interaction strength, from zero 
to a finite value. Such experiments have been performed 
with ultracold atoms, for example, Hung et oZ.^suddenly 
decreased the interaction strength in a weakly interacting 
atomic Bose-Einstein condensate . The resulting density 
oscillations were related to Sakharov oscillations in the 
early univers^^. 

An even closer connection can be made with the split¬ 
ting quench by Langen et al. in one-dimensional atomic 


condensateJi^. When a condensate is rapidly split in two 
parts, there is initially perfect phase coherence between 
them. However, at later times, the two parts start to 
develop a different phase. This dephasing due to inter¬ 
actions is entirely analogous to the one in our dynamical 
Casimir proposal, showing a light-cone-like emergence of 
thermal correlations. 

An important difference between polaritonic and 
atomic condensates concerns the ratio of the life time 
with respect to the characteristic time scale of the dy¬ 
namics. Whereas for ultracold atoms, this ratio is very 
large, in polariton systems losses are more important. 
Their theoretical modelling should therefore be carried 
out in an open system setting. This raises the interesting 
issue of the competition between losses and thermaliza- 
tion dynamics. 

We will treat the open system quantum dynamics 
within the truncated Wigner approximation, which is 
a po pular tool in both the study of conservative cold 
atom^i^^i^ as for lossy polariton system^. When the 
condensate depletion is small, the equations of motion 
can be linearized in the fluctuations, which is equivalent 
to the Bogoliubov approximation. 

In Sec. |nj we use this approximation allows to obtain 
analytical results for the first and second order coher¬ 
ence functions in the homogeneous system. For the in¬ 
homogeneous case instead, we perform in Sec. [TV] Monte 
Carlo simulations of the stochastic equations of motion. 
We show that for a large smooth pumping spot, local 
density approximation satisfactorily reproduces the first 
order coherence function. Conclusions are drawn in Sec. 
0 


II. THE MODEL 

When a microcavity is excited sufficiently close to the 
lower polariton branch and all the relevant energy scales 
(linewidth, interaction energy) are much smaller than the 
Rabi splitting, it is well justified to restrict the dynamics 
to the lower polariton branch. 

We consider a driven dissipative bosonic system, whose 
dynamics is governed by a master equation of the Lind- 
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blad type 


space 


(1) 

Here, the Hamiltonian H = Hp + Hp contains the free 
Bose gas dynamics of the polaritons 


Hp = J dx {x) 


Z^, 

2m 




il){x), (2) 


with m being the lower polariton effective mass and g the 
interaction strength, and it includes the external laser 
driving 


Hl= J dx [FLix,t)-ip^{x) + Fl{x,t)'>p{x)] 


( 3 ) 


(l){x,t) = cj)c{t) + ( 6 ) 

where L is the length of the one dimensional wire that 
we consider. The evolution of the condensate density is 
determined to be 

{(t>*c{t)(t>c{t)) = nc{t) = ndO) exp {--ftH) (7) 

and the equations of motion for the fluctuations 
are linearised in In terms of the vector 

^{k) = [4>{k), (j){—k)]'^ and the noise vector dE{k) = 
[dW{k), dW{—k)]'^, they read 


where Fp is the laser amplitude. 

The polariton losses, which depend on the linewidth 7 , 
are described by the dissipator I?(/o), that we take to be 
of Lindblad form 

'^{p) = J ^ [2i;{x)pi;'<{x) 

{x)'ip{x)p — pip^ {x)'ip{x)]. (4) 

For the quantum quench that we consider, we take 
the driving laser to be an ultra short pulse. When the 
pulse duration St is much shorter than all the other 
time scales of the dynamics, shortly after the pulse, 
the polariton held is in a coherent state with amplitude 
tpo = {ip{x,t = 0)) = j'^g^Fp{x,t)dt. The external laser 
drive then only sets the initial condition and does not 
affect the polariton dynamics, which is governed by the 
free Bose gas dynamics and the losses only. 

We will solve the master equation 0 within the trun¬ 
cated Wigner approximation (TWA), a method that is 
widely used for the simulation of weakly interacting one¬ 
dimensional atomic condensates. The addition of losses 
makes the TWA even a better approximation to the exact 
dynamics. 

The resulting stochastic equations of motion reacP^ 


ih d4>{x, t) 


O' 

+5l</>(a;.i)P (l){x,t)dt 


-k 


hrj 

4AV 


dW{x, t), 


( 5 ) 


where AH is the volume of a single cell of the discretized 
grid. Since the expectation values of the stochastic helds 
are equal to the symmetrised averages of the quantum 
fields, the TWA can be used to study the quantum fields. 


III. HOMOGENEOUS SYSTEM 

A. Bogoliubov approximation 

As long as the condensate depletion is small, the dy¬ 
namics can be treated in the linearised Bogoliubov ap¬ 
proximation. The field (p{x,t) is decomposed in Fourier 


ihd^{k) = B{k, t)^{k)dt -|- 



( 8 ) 


where the Bogoliubov matrix equals 


B{k, t) 


f e(fc) + gricit) - ^ gucit) 

V -gndt) -e(fc) - gndt) 



(9) 


and e(k) = fdk^/{2m). From the solution of these 
stochastic differential equations, we can compute the 
time evolution of the correlation functions. 


B. Momentum distribution and first order 
coherence 


The differential equation ^ for the stochastic helds 
can be solved exactly in the limit fc —> 0 , which gives for 
the momentum distribution 


( 10 ) 

For large momenta, we will resort to the sudden 
approximatiorfi^, which has yielded a good description 
of the average value of the momentum distribution: 


{iljHk,t)tp{k,t)) = 


gncjO)] 

hujB{k) 


sin [hhjB{k)t^ 




( 11 ) 


where hLOB{k) = \^€{k) [e(fc) -k 2gnc(0)] is the Bogoli¬ 
ubov dispersion. In the next section, we will calculate 
the second order coherence from the momentum distri¬ 
bution and the anomalous average {ipik, t)il}{—k, t)). The 
latter quantity is calculated following the same procedure 
as for the momentum distribution. Thus, we hrst deter¬ 
mine {'ip{k,t)'ipi—k,t)) for a system without decay, i.e. 
7 = 0 . and reintroduce the time dependence by letting 
the expectation values decay exponentially. This yields 




gnc(0)sinlhuiB(k)t] ^-yt/h n‘-i\ 


{[e(/c) -k gnc{0)] sin [fkjJB{k)t] + ihujB{k) cos [huJB{k)t\} . 
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The Fourier transform of the momentum distribution 
gives us the first order correlation in real space 




Vii’"’ (x, t)'ip{x, t)) (^/>t (a;', t)^ix', t)) 


(13) 


From this quantity, we can obtain the condensate frac¬ 
tion, the quantum depletion Snlric, 


Sn{t) ^ g^ne(O) 
Ticit) 72 


1 - ( ? + 1 
h 


(14) 


with 




! li _ 12! + li 

V g«c(o) [ \h ) 


- 1 - 1/2 


(15) 


and the coherence length ic 


4(0 = 2 . 1 /A:*(<), 


(16) 


where the factor 2.1 was determined numerically. 


C. Second order coherence in momentum space 

Since the particles are predicted to be produced in 
pairs with opposite momentum, we expect to find a cor¬ 
relation between polaritons with momentum k and those 
with momentum —k. Therefore, we will study the second 
order coherence in momentum space: 

g^'^\k,-k,t) = (17) 

(V:l'(fc,f)4l'(-fc,O'0(-fc,t)4(fc,O} 

(i/>t(fc, 00 (fc, t)){tp^-k, f) 4 (-fc, t)) ■ 

By applying Wick contraction, we can write the de¬ 
nominator as a product of quadratic expectation values. 
The nonzero terms are those containing the momentum 
distribution {ij}^ (k,t)'(p{k,t)) and the anomalous average 
{tlj^k,t)'4>\—k,t)). In terms of the stochastic fields, the 
expression becomes 

(i/>'i'(fc,t)4'i'(-fc,t)V’(-fc,t)4(fc,t)) = 

{(j)*{k, t)(j){k, t)) {(j)*{-k, t)(j){-k, t)) (18) 

+{(j)*{k,t)(j)*{-k,t)){(j){-k,t)(j){k,t)) 

-\{4'*{k,t)(j){k,t)) - ^{(l)*{-k,t)(l){-k,t)) + i. 

In the limit fc —>■ 0, an exact solution of equation (|^ can 
be found, which yields for the second order coherence 


lim —fc,t) = 2-|- (19) 

fe-fO 

(yf/h)^ exp {—2'yt/h) 

4 (gnc(0)/7)2 [(1 -|- 'yt/hY exp {—A'yt/h) — 2(1 -|- ^t/h) exp {—Z'yt/h) + exp (— 274 ^)]' 


The fraction in this expression diverges both at short 
times (t ^ ^/t)) when \\Ta.k^Qg^‘^\k,—k,t) k 2 + 
{gnc{0)t/h)~‘^ and for long times {t S> k/^), when 
limfc_>o g^^\k, —k,t) « 2 -I- {'j'^t/hgnc{0))‘^■ However, for 
good cavities with 7 <C gUciO), there is a large time win¬ 
dow when limfc_j.o —k, t) ~ 2. 

For large momenta, an expression for the second order 
coherence can be calculated from equation (181 and the 


results obtained with the sudden approximation ( 11 ) and 


(12). Combining these solutions and averaging out the 


oscillations afterwards, yields the following expression for 
the second order coherence 


g^'^\k,-k,t) =2 + 2 


hwsik) 

gnc{0)] 


( 20 ) 


Here, we recognise the inverse of the momentum distri¬ 
bution without the exponential decay. This analytical 
result can be well understood from the assumption that 
the particles are indeed produced in pairs. In this case, 
a polariton with momentum k will always be accompa¬ 
nied by a polariton with momentum —k. For large mo¬ 
menta, we expect very few particles. When the expected 


value is much smaller than one, only in few of the real¬ 
isations a polariton will be present and the cases with 
more than one polariton with the same momentum will 
be negligible. Therefore, we find that when the polaritons 
are produced in pairs with opposite momentum, when 
tp\k,t)'ijj{k,t) <C 1 , and neglecting the exponential de¬ 
cay, 

{tf}^k,t)tp\-k,t)'tp(-k,t)jp{k,t)) = {ip\k,t)tf}{k,t)). 

( 21 ) 

As a result, the second order coherence is the inverse of 
the momentum distribution 


g^^Hk,-k,t) 


large fc {tlj^k,t)tp{k,t))'' 


( 22 ) 


which is in close agreement with ( 20 ). 


For the homogeneous system, the numerical calcula¬ 
tions have been performed using the Green’s function 
method fronP. In terms of the Green’s function 


N 

Gkit,t') = Y[exip[-iAtBk{tj)/H], (23) 
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Figure 1. Second order coherence in momentum space for 
a homogeneous system calculated with the Green’s function 
method. Solid lines show the numerical results, whereas 
the dotted line represents eq. ( |20[ |. We have chosen g = 
0.01 ^mmeV, 7 = 0.05 meV, ft = 1, m = 1, gnc{0)/j = 10. 



Figure 2. Momentum distribution for a homogeneous system 
calculated with the Greens’ function method (blue solid) and 
Monte Carlo simulations (red dotted). We have chosen g = 
0.01 ^mmeV, 7 = 0.05 meV,!/ = 200/rm, ft = 1, m = 1, 
gnc{Q)h = 10 . 


the second order coherence can be written 


= 1 + 


Glit, 0)J ^ ^ + [Glit, 0)J ^ ^ [Gkit, 0)]i^2 






J 1,2 


(24) 


+ 


Gl{t,s)\^JGUt,s)],^^ 


X [{ip\k,t)'ip{k,t)){ip\-k,t)ilj{-k,t))] \ 


where the x, y of [Gfc(t, s)],^,, y indicate the matrix com¬ 
ponent. This results is depicted in figure together with 
the analytical expression eq. ( [^ derived from the sud¬ 
den approximation. It can be seen that the analytical ex¬ 
pression describes the overall behaviour of {k, —k,t) 
very well. The fast oscillations of the numerical results 
are expected to become more averaged in experimental 
data, which would then become closer to the analytically 
calculated average. 


IV. INHOMOGENEOUS SYTEM 
A. Monte Carlo simulation 

Although the Green’s function method provides a good 
description of the homogeneous system, it has some lim¬ 
itations. First, the interactions energy should not be too 
high, since the Bogoliubov approximation is no longer 
valid when the quantum depletion becomes too large. 
Secondly, the Green’s function method becomes cumber¬ 
some for inhomogeneous systems. In order to overcome 
these problems, we have also implemented a Monte Carlo 
simulation algorithrrfi^. 


In the truncated Wigner Monte Carlo algorithm, the 
expectation values are calculated by averaging over many 
realisations of the system. The results presented here are 
obtained from 10000 realisations. For the initial situa¬ 
tion, an average density is chosen and random noise is 
added to account for the stochastic nature of the fields. 
As opposed to the Green’s function method, where equa¬ 
tion (|^ was solved, the Monte Carlo algorithm used 
both the real and momentum space representation of the 
stochastic fields. The evolution due to interactions and 
decay, which are the time-dependent parts of the Hamil¬ 
tonian, is calculated in real space, whereas the effect of 
the kinetic term is calculated in momentum space. This 
method has the advantage that we do not have to dis¬ 
tinguish between condensate and excitations in the in¬ 
teraction term. In order for the algorithm to work, the 
time steps for which the evolution is calculated should 
be small, in order for the effect of sequentially calculating 
the evolution in real an momentum space to be small. For 
the homogeneous system we have used a system length 
of 200 fim and 0.4 ^m as the size of a unit cell. For the 
inhomogeneous systems, we adapted the length and grid 
size in order to have sufficient detail, the boundaries of 
the system distant enough with respect to the width of 
the Gaussian distribution, while keeping the number of 
grid points equal to 265 for computational efficiency. 

First, we verified the results of the Green’s function 
method for the homogeneous system, see Fig. Sec¬ 
ondly, we studied the effect of a Gaussian density distri¬ 
bution. For the second order coherence, the expectation 
value of four operators can be computed directly using 
the Monte Carlo simulations as 

{t{j\k,t)'tp^-k,t)'tlj{-k,t)-tp{k,t)) = 

{(t)*{k, t)4>* {-k, t)(j){-k, t)4>{k, t)) (25) 

{k,t)(l){k,t)) - i((^*(-fc,t)<))(-fc,t)) -b 
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Figure 3. Density calculated with Monte Carlo simula¬ 
tions, initial Gaussian profile given by exp(—with 
s = 100 ^m. We have chosen g = 0.01/rmmeV, 7 = 
0.05 meV, h — 1, m — 1, griij'Y — 10. 



x/s 


Figure 4. Density calculated with Monte Carlo simulations, 
initial Gaussian profile given by exp(—x^/s^), with s = 20 fim. 
We have chosen g = 0.01 /rmmeV, 7 = 0.05 meV, h = 1, m = 
1 , gm/^ = 10 . 


B. Numerical results 

The Monte Carlo simulation was applied to study sys¬ 
tems with an initial Gaussian density distribution, given 
by exp(—x^/s^). For the results presented here, the 
values s = 100 /rm and s = 20 ^m have been chosen and 
the initial central density rii = 50 ^m“^. For this density, 
the Bogoliubov approximation and truncated Wigner ap¬ 
proximation were still valid for the homogeneous sys¬ 
tem. For s = 100/xm, the behaviour in real space can 
be well understood. The density at the centre of the 
Gaussian decreases faster than the overall exponential 
decay, whereas the density at the sides of the distribu¬ 
tions shows a relative increase, see Fig. This would be 
expected from the repulsive interactions. As a result, the 
density distribution becomes more homogeneous. In the 
case of a smaller Gaussian distribution, where s = 20 /rm, 
this effect is even stronger. At the latest depicted time. 



Figure 5. First order spatial coherence g^^\x, —x,t) for an 
initial Gaussian distribution (exp(—x^/s^)) with width s = 
100 pm, calculated with Monte Carlo simulations. Ic is given 
in (|16[), and 5n n in (|14[). Same parameters as Fig. 


t = 39ps, the central region is very homogeneous, and 
shows a fast decay at the edges of the distribution. 

The first order coherence for the s = 100 pm system 
can be described by the results obtained from the homo¬ 
geneous systenP. The coherence length and the maximal 
depletion are close to that of the homogeneous system 
with the density that is equal to the maximal density 
of the Gaussian distribution rii. Since there is still a 
linear relation between the coherence length and the de¬ 
pletion, a simple correction to the numerical constants 
would give an even better description of these quanti¬ 
ties. The overall shape of {x, —x, t) is a direct result 
from the Gaussian shape of the density. For the homo¬ 
geneous system, a formula for the depletion at very large 
times as a function of the blueshift gndt) was derived: 
Sn/n « 0.77 5 /(^ 7 ). In the inhomogeneous system the 
healing length ^ = h/^Jmgricifi) becomes position depen¬ 
dent. Consequently, the final depletion will also depend 
on the position. When nc(0) is simply replaced by the 
initial Gaussian distribution, exp(—x^/s^), the black 
dashed lines from Figs. and are found. For the wider 
Gaussian, with s = 100 pm, this describes the behaviour 
of the first order coherence very well. Therefore, the phe¬ 
nomenon that the coherence goes back to one for large 
distances is due to the small density which leads to a 
smaller depletion. However, for the smaller Gaussian, 
where s = 20 pm, the local density approximation is no 
longer valid. Nevertheless, the initial decay of the first 
order coherence, until x « 20 pm is still described by the 
linear relation found in the homogeneous case. 

In momentum space, the system with a Gaussian den¬ 
sity is very different from a homogeneous system. This 
is in accordance with expectations, since for an initial 
Gaussian distribution, many momentum states are occu¬ 
pied from the start. Moreover, the inhomogeneous den¬ 
sity profile leads to an expulsion of the polaritons away 
from the region where the were created. This acceler- 
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Figure 6. First order spatial coherence g^^\x,—x,t) for 
an initial Gaussian distribution (exp(—x^/s^)) with width 
s = 20 pm, calculated with Monte Carlo simulations. Same 
parameters as Fig. 



Figure 7. Momentum distribution for an initial Gaussian dis¬ 
tribution (exp(—with width s = 100 pm, calculated 
with Monte Carlo simulations. Same parameters as Fig. 

ation corresponds to a shift of the momentum distribu¬ 
tion. Hence, momentum conservation in the interactions 
no longer results in pair production of polaritons with 
opposite momentum. In Figs. an di we see an increase 
of the particle number with respect to the homogeneous 
case at small yet finite momentum, which indicates that 
the distribution is expanding, which is indeed what was 
seen in the evolution of the density. The wider and the 
smaller distribution both display very similar behavior, 
where the smaller density distribution has a peak in the 
momentum distribution at larger momenta, as compared 
to the wider density distribution. At larger momenta, the 
momentum distribution follows that of the homogeneous 
system more closely for the s = 100 pm case. Neverthe¬ 
less, the particle number at these momenta is very small. 

The second order coherence, depicted in Figs. and 
|10[ is even more different from the homogeneous case. 
For small momenta, the second order coherence is close 



Figure 8. Momentum distribution for an initial Gaussian dis¬ 
tribution (exp(—with width s = 20pm, calculated 
with Monte Carlo simulations. Same parameters as Fig. 



Figure 9. Second order coherence in momentum space for 
initial Gaussian distribution (exp(—with width s = 
100 pm, calculated with Monte Carlo simulations. Same pa¬ 
rameters as Fig. 

to one, which is the value corresponds to a coherent sys¬ 
tem. We see that the second order coherence remains 
close to one, even for momenta at which the momentum 
distribution seems to resemble the homogeneous case for 
the wider density distribution. This suggests that even 
for these momenta, most particles are part of the conden¬ 
sate. For the system with s = 20 pm, the second order 
coherence, is either very close to one, or it contains too 
much noise for a good description. 


V. CONCLUSIONS 

We have studied a quantum quench consisting of a 
sudden injection of polaritons in a microcavity. Both 
a homogeneous and a Gaussian initial density distribu¬ 
tion have been examined. The homogeneous case has 
been related to the dynamical Casimir effect previously. 
Where the correlation functions in the homogeneous case 
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could still be calculated analytically, we performed trun¬ 
cated Wigner Monte Carlo simulations for the Gaussian 
excitation pulse. The first order spatial coherence is well 
approximated by the local density approximation for a 
sufficiently wide pulse. The second order coherence in 
momentum space evidences the production of excitations 
in pairs. For a system with an initial Gaussian density 
distribution, multiple momentum states are significantly 
occupied from the start. Here, the second order coher¬ 
ence indicates that many of these particles are coherent, 
so that evidence of quantum correlations is highly sup¬ 
pressed. 


Figure 10. Second order coherence in momentum space 
for initial Gaussian distribution (exp(—x^/s^)) with width 
s = 20 pm, calculated with Monte Carlo simulations. Same 
parameters as Fig. 
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